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Abstract 

Background: Huanglongbing (HLB) is arguably the most destructive disease for the citrus industry. HLB is caused 
by infection of the bacterium, Candidatus Liberibacter spp. Several citrus GeneChip studies have revealed thousands 
of genes that are up- or down-regulated by infection with Ca. Liberibacter asiaticus. However, whether and how 
these host genes act to protect against HLB remains poorly understood. 

Results: As a first step towards a mechanistic view of citrus in response to the HLB bacterial infection, we 
performed a comparative transcriptome analysis and found that a total of 21 Probesets are commonly up-regulated 
by the HLB bacterial infection. In addition, a number of genes are likely regulated specifically at early, late or very 
late stages of the infection. Furthermore, using Pearson correlation coefficient-based gene coexpression analysis, we 
constructed a citrus HLB response network consisting of 3,507 Probesets and 56,287 interactions. Genes involved in 
carbohydrate and nitrogen metabolic processes, transport, defense, signaling and hormone response were 
overrepresented in the HLB response network and the subnetworks for these processes were constructed. Analysis 
of the defense and hormone response subnetworks indicates that hormone response is interconnected with 
defense response. In addition, mapping the commonly up-regulated HLB responsive genes into the HLB response 
network resulted in a core subnetwork where transport plays a key role in the citrus response to the HLB bacterial 
infection. Moreover, analysis of a phloem protein subnetwork indicates a role for this protein and zinc transporters 
or zinc-binding proteins in the citrus HLB defense response. 

Conclusion: Through integrating transcriptome comparison and gene coexpression network analysis, we have 
provided for the first time a systems view of citrus in response to the Ca. Liberibacter spp. infection causing HLB. 



Background 

Orange juice is among the largest beverage industries in 
the world. Sweet orange is mainly produced in the sub- 
tropical areas in the countries of China, US and Brazil and 
the Mediterranean basin regions. Sweet orange (Citrus 
sinensis) belongs to the Citrus genus that includes several 
other species such as tangerine, mandarin and grapefruit. 
In horticultural practice, citrus is asexually propagated 
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through grafting the scion onto the stock which is grown 
through the seeds. The scion has been bred for the desired 
traits of fruit quality while the stock is mostly selected for 
supporting the optimal scion growth and increased resist- 
ance to biotic and abiotic stresses. 

Among the major biotic factors which frequently chal- 
lenge tree growth and fruit development, Huanglongbing 
(HLB) or called citrus greening is one of the most de- 
structive diseases. HLB was first reported in 1919 in 
southern China, and very recently it has been reported 
in almost all major citrus production areas [1-3]. For ex- 
ample, in Florida alone, HLB has caused the loss of se- 
veral billion dollars since 2005 when HLB was first 
reported, ranging from 30-100% of loss in fruit production 
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in citrus groves [3]. HLB is caused by infection of the bac- 
terium, Candidatus Liberibacter spp., which is spread to 
plants via the vector Asian citrus psyllid (Diaphorina citri) 
or through grafting of a diseased shoot. The HLB bacter- 
ium has three species, Ca. Liberibacter asiaticus (Las), Ca. 
Liberibacter africanus (Laf) and Ca. Liberibacter ameri- 
canus (Lam). The genome of the Las species was recently 
published, with a size of approximately 1.23 Mb [4]. It has 
been generally accepted that, after infection or inoculation, 
the HLB bacteria migrate through phloem and, by accu- 
mulating there, causes the formation of sieve plug [2,3,5-7]. 
Consequently, the transport of nutrients (such as carbohy- 
drates and amino acids) from the source leaves to various 
sinks (such as young leaves, fruits and roots) are compro- 
mised or even blocked in severely infected plants, leading 
to the alterations in carbohydrate metabolism for meta- 
bolic flow and exhibiting such phenotypes as yellow and 
blotchy mottles on leaves, variegated fruits and poor root 
growth [6-8]. Because of the huge impact of HLB in the 
citrus industry, plant pathologists and horticulturists have 
long sought after the HLB resistance mechanism in citrus. 

A recent survey suggests the existence of genetic varia- 
tions among different citrus species, varieties and stocks 
[9]. In general, mandarin, sweet orange and grapefruit are 
relatively more susceptible to the HLB bacterial infection, 
while sour orange, lemon, lime, and citrange (a hybrid of 
sweet orange and the stock C. trifoliata) are less suscep- 
tible. This raises the possibility that HLB resistance can be 
achieved through genetic means. Nevertheless, breeding 
for the HLB resistance through crossing will be a daunting 
task, given the complex genetic backgrounds, the nature 
of asexual propagation and the relatively long juvenile 
period for citrus. Therefore, many researchers have 
turned their attentions to finding the target genes that 
are required or critical for the citrus host response to the 
HLB bacteria. Transcriptome analysis has been used as a 
straight forward approach to identify the genes whose ex- 
pression is altered in citrus leaves in response to the 
HLB inoculation [5,10-12]. These studies led to the iden- 
tification of several hundred or thousand genes that are 
up- or down-regulated by the HLB bacterial infection. 
The majority of these genes can be grouped into 



metabolism, transport and response to stimulus. However, 
these studies varied significantly in terms of study design 
(including experimental materials and infection stages) and 
data analysis (for example, different statistical methods and 
various fold-change cutoffs). Furthermore, there is a lack of 
comparison of the results from these different experiments. 
In addition, how these HLB bacterium-regulated genes are 
connected in a system remains unknown. 

To provide a systems view of citrus response to the HLB 
bacterial infection, we first performed a comparative study 
of the previously reported transcriptome datasets. Our 
results show that there are 21 probe sets (representing up 
to 19 genes) are commonly up-regulated and a number of 
genes that are specific to early, late or very late stages of in- 
oculation. Furthermore, using the Pearson correlation coef- 
ficient (Pec) -based unweighted gene coexpression analysis, 
we constructed an HLB response network. This citrus gene 
coexpression network consists of 3,507 Probesets and 
56,857 interactions. We then mapped certain categories of 
the HLB responsive genes to the HLB response network, 
resulting in the formation of several important subnet- 
works including metabolism, transport, signaling, defense 
response and hormone response. Taken together, through 
comparative transcriptome analysis and construction of a 
citrus gene coexpression analysis, we have provided a sys- 
tems view of citrus response to the Ca. Liberibacter infec- 
tion causing HLB. 

Results 

An overview of comparative analysis of HLB 
transcriptomes 

To perform a comparative transcriptome study, we 
decided to use the same data pre-processing and statis- 
tical analysis methods and the same selection criteria for 
the identification of HLB significantly regulated genes. 
Two sets of the citrus Affymetrix GeneChip data derived 
from very recent publications [11,12] were retrieved 
from the NCBI Gene Expression Omnibus (GEO) data- 
base, while the data for the two earlier reports were pro- 
vided by Drs. Bowman [10] and Wang [5], respectively. 
These four reports represent six different studies that 
can be used for individual comparisons, with a total of 



Table 1 A list of citrus GeneChip studies used in this analysis 



Studies 


Varieties 


Weeks after 
inoculation 


Ref. 


Stages 


No. of significantly regulated Probesets reported in this study 
Up-regulated Down-regulated 


1 


Valencia 


5-9 


[10] 


Early 


158 


62 


2 


Valencia 


13-17 


[10] 


Late 


394 


285 


3 


Sweet orange 


35 


[5] 


Very late 


942 


1060 


4 


Madam Vinous 


30 


[11] 


Very late 


1955 


2125 


5 


Cleopatra 


30-32 


[12] 


Very late 


1110 


626 


6 


US-897 


30-32 


[12] 


Very late 


261 


39 
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34 arrays (see Table 1). In these studies, genome-wide 
gene expression was profiled from the citrus leaves 
inoculated by the HLB bacterium Las. However, these 
six studies can be categorized into three distinct HLB dis- 
ease stages (see Table 1). Because the three studies used 
the leaf samples 30-35 weeks after inoculation, we arbi- 
trarily called this "very late" stage by following the defin- 
ition of "early" (5-9 weeks) and "late" stages (13-17 weeks) 
described in the first HLB transcriptome study [10]. 

The citrus GeneChip contains a total of 30,173 Probe- 
sets. Because the Affymetrix company has not provided a 
comprehensive annotation analysis for those Probesets, it 
is not known how many unique citrus genes are actually 
represented in the chip. Therefore, we decided to analyze 
the number of Probesets that are significantly regulated 
in response to HLB. The data pre-processing was 
described in Methods. In brief, those Probesets with the 
calls of present (P) or marginal (M) in at least two chips 
in each of the four reports were included in our analysis. 
For the identification of significantly regulated genes, the 
adjusted LPE approach was used because of its power in 
analyzing small samples [13]. In our analysis, a two-fold 
cutoff was used, resulting in various numbers of genes 
that were either up- or down-regulated in each of the six 
studies (see Table 1). The HLB regulated genes for each 
study were listed in Additional file 1. If the genes signifi- 
cantly regulated in at least one study were added to- 
gether, we found that a total of 3,345 and 3,230 Probesets 
were up-regulated and down-regulated, respectively. 
These Probesets are called "HLB responsive genes" in 
this study (Additional file 1). The percentage of "HLB re- 
sponsive genes" identified in this comparative analysis 
(22%) is similar to that of the bacterial pathogen respon- 
sive genes in Arabidopsis [14]. This indicates that either 
the disease response mechanism could be somehow con- 
served or these four reports have probably identified 
most of the HLB responsive genes in the citrus genome. 
Surprisingly, the study-wise comparison showed that the 
proportion of the significantly regulated genes over- 
lapped in two of the six studies varied dramatically (from 
0.05% to 75%) (see Additional file 2). 

Comparative studies reveal commonly regulated and 
stage-specifically regulated genes by HLB 

Despite our finding that only a small proportion of Pro- 
besets are significantly regulated in any of two studies, 
we reasoned that those Probesets (and their represented 
genes) commonly regulated in all of the studies may rep- 
resent either a common core pathway or default pathway 
in response to the Las infection. We first found a total 
of 13 Probesets that are commonly up-regulated in all of 
the six studies (see Table 2), representing only 0.4% (14 
out of 3,426) of the HLB up-regulated genes. However, 
the number of Probesets significantly regulated in any of 



five studies increased to 42. It is possible that in the ab- 
sence of the HLB bacterial challenge some of the HLB 
up-regulated genes already had higher transcript levels 
in the relatively resistant germplasm US-897 compared 
to the relatively susceptible mandarin Cleopatra and thus 
they were not up-regulated any more in US-897 in re- 
sponse to the Las infection; however, they could be sig- 
nificantly regulated in all other four studies. We did 
identify a total of eight Probesets for this type of expres- 
sion pattern and consequently they were also added to 
the list of the HLB commonly regulated genes (see 
Table 2). Surprisingly, there was no Probeset commonly 
down-regulated in all of the six studies and only one 
Probeset that is significantly down-regulated in five stud- 
ies. This Probeset, Cit.l8719.1.Sl_at, is annotated to en- 
code a gene similar to Arabidopsis AT5G18600-encoded 
glutaredoxin family protein involved in cell redox 
homeostasis. Gene Ontology (GO) analysis of the subset 
of 21 commonly up-regulated Probesets indicates that 
metabolism (5 Probesets), transport (5 Probesets), hor- 
mone responses (3 Probesets) and unknown processes 
(7 Probesets) are the largest groups. The three Probe- 
sets representing the genes involved in hormone re- 
sponse indicate that gibberelic acid (GA), abscisic acid 
(ABA), auxin, ethylene and jasmonic acid (JA) may 
have certain role in mediating the citrus response to 
HLB. Interestingly, three Probesets (Cit.30071.1.Sl_at, 
Cit.33707.1.Sl_at, and Cit.3534.1.Sl_s_at) belonging to 
the category of "unknown process" might also be 
involved in ethylene response as they exhibit the high- 
est homology to genes that are associated with ethylene 
response using the manual BLAST search. In addition, 
there is a transcription factor gene represented by the 
Probeset, Cit.l2214.1.Sl_s_at (an NAC domain contain- 
ing protein called NAC096 in Arabidopsis) and another 
putative RAP2.4-like ethylene transcription factor repre- 
sented by Cit.3534.1.Sl_s_at. Taken together, the exist- 
ence of these commonly up-regulated genes strongly 
indicates that metabolism, transport, hormone response 
and transcriptional regulation play a critical role and 
may define the default or basal pathways in citrus dur- 
ing the whole process of the Las infection. 

In contrast to the commonly regulated genes in HLB 
response, we found various numbers of stage-specifically 
regulated genes as this group of genes were only regu- 
lated at a particular stage. There are 27 and 7 Probesets 
that are respectively up-regulated and down-regulated 
only at early stage, and thus they belong to the group of 
early stage-specific genes (Table 3). In this early stage- 
specific group, the genes encoding components involved 
in cell wall metabolism and transcription are of particular 
interest. First, there are six Probesets that could poten- 
tially represent the genes involved in cell wall biogenesis 
or property (Cit.l7310.1.Sl_s_at, Cit.l7450.1.Sl_s_at, 
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Table 2 A list of commonly up-regulated HLB responsive genes in various studies 



Probeset ID AtGID Arabidopsis gene annotation GO term 



LJ L. I UUjZ. 


1 1 
1 .J 1 


_X_3t 


AT1 C7^7^ri 


GA-responsive GAST1 homolog; BR, GA and ABA responsive expression; 
located in cell wall; GASA1 (GAST1 PROTEIN H0M0L0G 1) 


Response to stimulus 


Cit.12214. 


1.S1 


_s_at 


AT5G46590 


NAC096 (Arabidopsis NAC domain containing protein 96); transcription factor 


Transcription 


Cit.13437. 


1.S1 


_s_at 


AT4G39210 


APL3 (sugar-inducible); glucose-1- phosphate adenylyltransferase 


Carbohydrate metabolic process 


Cit.17561. 


1.S1 


_s_at 


AT5G05365 


Metal ion binding 


Transport 


r\t 1 Qf^7/i 

LJ L. I ~>0/ H. 


1 ^ 1 
I .J I 


_s_at 


A I I ^HJ 


L0X2 (LIPOXYGENASE 2); lipoxygenase; response to JA, bacterium, fungus 
and wounding; JA biosynthetic process 


Response to stimulus 


Cit.25840. 


1.S1 


_s_at 


AT1G 17840 


WBC1 1 (WHITE-BROWN COMPLEX H0M0L0G PROTEIN 11); ATPase, coupled 
to transmembrane movement of substances / fatty acid transporter; 

rccmnco ~\~/~\ ctmcc nr>r\ ARA 

response lu suess anu ada 


Transport; Response to stimulus 


LIToUU/ 1 . 


1 C 1 

I .J I 


_3t 


A I jIjUz joU 


Unknown protein (auxin-ethylene interaction, gemivirus induced) 


Unknown process 


LIT.jj/U/ . 


1 CI 

I .J I 


_at 


AT/iri7onri 
A 1 4oz/yuu 


CCT motif family protein (ethylene, gemivirus-induced) 


Unknown process 


Lit. /Odd. I 


C 1 
o l_ 


_at 


A I D\3 I jodU 


Phospholipase/carboxylesterase fami y protein 


Other metabolic process 


Cit. 5367.1 


C 1 


.at 




Transcript Assignment: DN622378; Human protein tyrosine phosphatase 
type IVA 


Unknown process 


Ctf ^^/i 1 
L.jjj't. I 


J I _ 


_s_at 




Grape RAP2-4- ike ethylene transcription factor 


Unknown process 


Cit.29171. 


1.S1 


_at 




UCRPT02_65A04_b Poncirus trifoliate! Roots with Iron Deficiency cDNA clone 


Unknown process 


Cit.19313. 


1.S1 


_at 




Transcript Assignment: CX308038; Hypothetical protein 


Unknown process 


Cit.11459. 


1.S1 


_s_at 


AT1G05300 


ZIP5; cation/metal ion transmembrane transporter 


Transport 


Cit. 11460. 


1.S1 


_at 


AT1G05300 


ZIP5; cation/metal ion transmembrane transporter 


Transport 


Cit.1 21 72. 


1.S1 


_at 


AT4G35160 


O-methyltransferase family 2 protein 


Other metabolic process 


Cit.15630. 


1.S1 


_at 


AT4G33040 


Glutaredoxin family protein 


Oxidation-reduction process 


Cit.1 7155. 


1.S1 


_at 


AT1G02850 


BGLU1 1 (BETA GLUCOSIDASE 11); catalytic/ cation binding / hydrolase, 
hydrolyzing 0-glycosyl compounds 


Carbohydrate metabolic process 


Cit.1 7155. 


1.S1 


_s_at 


AT1G02850 


BGLU1 1 (BETA GLUCOSIDASE 11); catalytic/ cation binding / hydrolase, 
hydrolyzing 0-glycosyl compounds 


Carbohydrate metabolic process 


Cit.3171.1 


si. 


x_at 


AT1G 18880 


Proton-dependent oligopeptide transport (POT) family protein 


Transport 


Cit.23598. 


1.S1 


_s_at 






Unknown process 



The first section contains 13 ProbeSets up-regulated in all six studies. The second section includes eight Probesets up-regulated in five studies (except Study 6, 
see Table 1) but their expression was higher in US-897 than in Cleopatra. AtGID, Arabidopsis gene ID. 



Cit.23171.1.Sl_x_at, Cit.33335.1.Sl_s_at, Cit.35345.1. 
Sl_s_at, and Cit.30425.1.Sl_x_at). Second, five Probe- 
sets represent the transcription factors homologous to 
Arabidopsis ERF5 (Ethylene responsive element-binding 
factor), ATAF1 and ARR9 (a two-component response 
regulator). In addition, Cit.l6537.1.Sl_at represents a 
GCN5-related N-acetyltransferase family protein, which 
might be involved in global transcriptional control 
through chromatin remodeling. This result implies that 
transcriptional control and cell wall property regulation 
are among the early events in citrus in response to the 
HLB bacterial attack. 

In addition, 103 up-regulated and 74 down-regulated 
Probesets are specific to the "late" stage of Las infection 
(Additional file 3). Interestingly, these Probesets repre- 
sent some genes that belong to the categories of metabol- 
ism of carbohydrate, nitrogen and lipids, hormone IAA 
metabolism, response to chemical stimulus, endomem- 
brane systems and extracellular regions. In addition, 



while several genes involved in cell wall property 
regulation (such as expansin) are up-regulated, some 
genes encoding transcription factors (such as MYB15/ 
52/103) and protein kinases are down-regulated. The 
most striking feature is that only seven Probesets rep- 
resent the very late stage-specific genes (Additional 
file 3). These include the genes that are most closely 
related to Arabidopsis C domain containing protein 
71 (a transcription factor), a copper-binding family 
protein, a trypsin and protease inhibitor family pro- 
tein/Kunitz family protein, a myosin heavy chain- 
related protein, two basic chitinase (CHIB) and one 
unknown protein encoded by Atlg42430. The small 
number of genes belonging to this very late stage- 
specific category is likely due to the various experi- 
mental conditions because only 26 Probesets are 
commonly up- or down-regulated even in the four 
studies within the same very late stage of Las infec- 
tion (Additional file 1). Nevertheless, as this group of 
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Table 3 A list of early stage-specific HLB responsive genes 



Probeset ID 




FC 


AtGID 


Arabidopsis gene annotation 


Up-regulated 








Cit.1 2388.1. S1 


_at 


2.1 


AT2G40095 


unknown protein 


Cit.1 3969.1. S1 


_s_at 


2.1 


AT5G52640 


ATHSP90.1 (HEAT SHOCK PROTEIN 90.1); ATP binding / unfolded protein binding 


Cit.1 5958.1 .SI 


_at 


2.3 


AT2G26560 


PLA2A (PHOSPHOLIPASE A 2A); lipase/ nutrient reservoir 


Cit.16156.1.S1 


_at 


2.7 


AT5G41330 


potassium channel tetramerisation domain-containing protein 


Cit.16537.1.51 


_at 


2.4 


AT2G32030 


GCN5-related N-acetyltransferase (GT) family protein 


Cit.1 71 87.1. S1 


_at 


2.7 


AT1G05575 


unknown protein 


Cit.1 731 0.1 .SI 


_s_at 


2.2 


AT1G32170 


XTR4 (XYLOGLUCAN ENDOTRANSGLYCOSYLASE 4); hydrolase, acting on glycosyl bonds 


Tit 1 74.SD 1 ^1 


s at 


2.7 


AT2G46330 


AGP 16 (ARABINOGALACTAN PROTEIN 16) 


LJ L. 1 O 1 DO. 1 .3 


_at 


2.2 


AT3G02840 


Immediate-early fungal elicitor family protein 


r\t 1 Q&n?. 1 c. 1 

LJ L. I oOUj. 1 .3 


_at 


2.5 






LJ L. I oOUj. 1 .3 


_s_at 


2.1 






fit 1 Qf.9, 1 ^1 

\~\ 1. 1 yoo. i .3 i _ 


_s at 


2.4 


AT1G61680 


TPS 14 (TERPENE SYNTHASE 14); S-lilool synthase 


C\t 711^ 
LI L.Z3 I / 1 . 1 .3 


x at 


2.0 


AT3G22120 


CWLP (CELL WALL-PLASMA MEMBRANE LINKER PROTEIN); lipid binding 


fit- 1/1Q7Q 1 C 

LiLz^y/y. i .3 


_at 


2.5 






LI L.ZyDU'H-. 1 .3 


_s_at 


3.8 


AT4G20820 


FAD-binding domain-containing protein 


Cit.29533.1.S1 


_s_at 


2.6 


AT5G47230 


ERF5 (ETHYLENE RESPONSIVE ELEMENT BINDING FACTOR 5); D binding / transcription activator/ 
transcription factor 


LlL.3 I yO'f. 1 .3 


_at 


2.2 






Tit 1 ^ 
LI L.33Z03. ! .3 


_at 


2.3 






Lit. DDDDD. \ .3 


_s_at 


2.4 


AT3G22120 


CWLP (CELL WALL-PLASMA MEMBRANE LINKER PROTEIN); lipid binding 


Cit.34429.1.S1 


_s_at 


2.8 


AT5G08350 


GRAM domain-containing protein / ABA-responsive protein-related 


Cit.3449.1.S1_ 


.at 


2.3 


AT1G01720 


ATAF1; transcription activator/ transcription factor 


Cit.3449.1 .S1_ 


.s_at 


2.4 


AT1G01720 


ATAF1; transcription activator/ transcription factor 


Cit.3450.1 .S1_ 


at 


2.5 


AT1G01720 


ATAF1; transcription activator/ transcription factor 


Cit.35345.1.S1 


_s_at 


2.1 


AT3G28340 


GATL10 (Galacturonosyltransferase-like 10); polygalacturote 4-alpha-galacturonosyltransferase/ 
transferase, transferring hexo 


Cit.4259.1 .S1_ 


at 


2.2 


AT1G 14870 


Response to oxidative stress; LOCATED IN: plasma membrane 


Cit.4690.1 .S1_ 


at 


2.8 


AT5G07990 


TT7 (TRANSPARENT TESTA 7); fiavonoid 3-monooxygese/ oxygen binding 


Cit.9944.1 .S1_ 


x_at 


2.1 


AT4G37990 


ELI3-2 (ELIGTOR-ACTIVATED GENE 3-2); aryl-alcohol dehydrogese/ mannitol dehydrogese 


Down-regulated 








Cit.1 6792.1. S1 


_at 


2.1 


AT3G57040 


ARR9 (RESPONSE REGULATOR 9); transcription regulator/two-component response regulator 


Cit.1 8730.1. S1 


_at 


2.1 






Cit.23498.1.S1 


_s_at 


2.2 


AT3G 12750 


ZIP1 (ZINC TRANSPORTER 1 PRECURSOR); zinc ion transmembrane transporter 


Cit.30425.1.S1 


_x_at 


2.2 


AT3G 12500 


AthCHIB (Arabidopsis BASIC CHITINASE); chitinase 


Cit.38866.1.S1 


_at 


2.5 






Cit.5121.1.S1_ 


.at 


2.3 


AT5G57685 


AtGDU3 (Arabidopsis GLUTAMINE DUMPER 3) 


Crt.7043.1i1. 


at 


3.3 


AT1G 17860 


Trypsin and protease inhibitor family protein/Kunitz family protein 



FC, fold changes; AtGID, Arabidopsis gene ID. 



genes were identified from four studies specifically at 
the very late stage compared to only one study for 
early and late stages respectively, they could be more 
reliable than groups of early- and late-stage-specific 
genes. 



Construction and characterization of gene coexpression 
network for citrus response to HLB 

To provide a systems view of citrus host response to the 
HLB bacterial infection, the Pearson correlation coeffi- 
cient (Pec) method was used to infer the gene 
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coexpression network using the four datasets reported in 
the four transcriptomic studies [5,10-12]. A total of 
10,668 Probesets, which are present (with P or M 
calls) in at least two chips of the transcriptomic stud- 
ies with strong expression and/or belong to the group 
of the HLB responsive genes (Additional file 1), were 
used for network analysis (see Methods). This number 
represents 35% of 30,173 Probesets in the citrus Gene- 
Chip. Pec was computed between each pair of these 
Probesets. A Pec threshold of 0.93 was selected, based 
on the overall consideration of changes in node num- 
ber, edge number and network density (Additional file 
4). The signed scale free R 2 plot analysis (Additional 
file 4) suggests that this selection has a good scale-free 
topology fit, as the R 2 value of 0.85 indicates that the 
topology of the HLB response network is quite similar 
to most biological networks. The resulting citrus gene 
coexpression network (called the "HLB response net- 
work" in this report) contains 3,507 nodes (Probesets) 
with 56,287 edges (interactions or connections) (see 
Figure 1 and Additional file 5). 



We next determined the robustness of our network 
across each dataset using the cross-validation approach. 
We randomly left out one dataset and reconstructed the 
gene co-expression networks using the remaining three 
datasets. The resulting four networks were then compared 
to the network based on all four datasets in terms of net- 
work connectivity rank of each gene according to the sug- 
gestion described elsewhere [15]. There were strong, highly 
significant connectivity correlations (R = 0.75, 0.83, 0.83, 
and 0.80 respectively, each associated with a p-value of less 
than 2.2e-16) between the network based on all four data- 
sets and the ones reconstructed from any combination of 
the three datasets. This suggests a high degree of preserva- 
tion of gene co-expression patterns across the networks 
based on different datasets. 

We then analyzed in detail the characteristics of the 
HLB response network. First, the frequency distribution of 
edges for each node was determined. As shown in Figure 2, 
the network contains 860 Probesets that are orphan nodes, 
400 Probesets that have only one interaction, and the ma- 
jority of the nodes (2,247 Probesets or 64%) that have at 





o 

o 



. Niirogmft amino »c*d 

" iii.i..-. 

Transport 
Otltntt rtioom* 
Setting 

Hormone rpspome 



Figure 1 HLB response network. The citrus HLB response network is displayed using Cytoscape. A white circle represents a node (Probeset), 
and a black edge connecting the two nodes indicates an interaction (coexpression) with a Pec value of > =0.93. Six categories of major biological 
processes are color-coded. 
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Figure 2 Summary of the network characteristics. The 

distribution of number of nodes (Probesets) in different ranges 
of edge numbers (1; 2; 3-10; 11-50, 51-100; 101-150; 151-200; 
201-250; 251-300; and 300-369) in the HLB response network is 
presented. 



least three interactions (see Figure 2) and, by following 
Geisler-Lee et al. [16], are called hubs in this paper. 
Among the 2,247 Probesets (or called hubs), the majority 
have 3-100 edges (1,880 Probesets or 84%), and the 
remaining 345 Probesets (15%) have 101-300 interactions, 
while only 1% (22 Probesets) have more than 300 interac- 
tions. Overall, the mean number of interactions for each 
Probeset is 16, with the maximum of interactions being 369 
(Cit.4987.1.Sl_s_at). Cit.4987.1.Sl_s_at represents a gene 
most closely related to Arabidopsis SYP71 encoding a plant 
syntaxin which functions as a plasma membrane- associated 
protein transporter. 

Second, we performed a GO enrichment analysis for 
the Probesets in the HLB response network. Among 
30,173 Probesets, 22,775 (or 75%) have the Arabidopsis 
gene ID (AtGID) as their closest orthologs or homologs 
(Additional file 6). Therefore, these Probesets were 
assigned GO terms based on the most recent Arabidopsis 
GO assignment. The remaining Probesets were given three 
general GO terms: "biological process" (GO:0008150), 
"molecular function" (GO:0003674), and "cellular compo- 
nent" (GO:0005575). GO enrichment analysis using the 
hypergeometric statistical method with the Hochberg false 
discovery rate (FDR) adjustment showed that many GO 
terms were overrespresented in the HLB response network. 
Among the overrepresented GO terms (Additional file 5), 
the nodes belonging to the following six categories were 
color coded in the HLB response network (see Figure 1): 
carbohydrate metabolic process (213 Probesets), nitrogen 
and amino acid metabolic process (207), transport (311), 
defense response (175), hormone response (238) and sig- 
naling (218). The nodes for each of these six categories, to- 
gether with the nodes belonging to some highly 
overrepresented GO terms such as response to stress, lipid 
metabolic process, cell wall and membrane part, were 



listed in Additional file 7. The p-values of the overrepre- 
sented GO terms were listed in Additional file 5. 

We also performed a GO enrichment analysis for 
the hub genes. We arbitrarily divided the 2,247 hubs 
into two categories: minor hubs (3-99 interactions) 
and major hubs (> = 100 interactions) and their overre- 
presented GO terms were summarized in Additional 
file 8. The major hubs have 13 overrepresented GO 
terms: carbohydrate metabolic process, primary meta- 
bolic process, metabolic process, secondary metabolic 
process, lipid metabolic process, cellular amino acid and 
derivative metabolic process, cellular process, localization, 
transport, establishment of localization, regulation of ana- 
tomical structure size, regulation of cell size, and regulation 
of cellular component size. In addition to these 13 GO 
terms, the minor hubs have 16 additional overrepresented 
GO terms, such as response to stimulus (including biotic 
and abiotic), response to stress, regulation of biological 
quality and signal transduction. 

Analysis of the defense and hormone response 
subnetworks 

Given the importance of carbon and nitrogen metabol- 
ism, transport, signaling, defense response and hormone 
response in the citrus response to the HLB bacterial in- 
fection and in general plant defense response, the sub- 
networks for these six categories were constructed by 
mapping the Probesets belonging to these categories into 
the HLB response network. The resulting edges (interac- 
tions) were listed in Additional file 7. 

We first analyzed the HLB defense subnetwork. As 
shown in Figure 3A, the Probesets representing defense, 
hormone response and signaling were color-coded. Clearly, 
the large hubs belong to the categories of defense and hor- 
monal responses but not signaling. Interestingly, several of 
the hormone response hubs are also defense response hubs 
as these hubs are involved in both responses (Figure 3B). 
For example, the Probesets Cit.ll529.1.Sl_s_at and 
Cit.ll530.1.Sl_at represent a transcription factor closest to 
the Arabidopsis At2G37630-encoded AS1 (Asymmetric 
leaves 1), which is annotated as both response to fungus, 
virus, bacterium and salt stresses and response to hor- 
mones such as auxin, GA, SA and JA [17-19]. Interestingly, 
these hubs were connected to other large defense hubs 
such as Cit.ll94.1.Sl_s_at, which represents a lipid trans- 
porter closest to Arabidopsis DIR1 (Defective induced re- 
sistance 1), Cit.3826.1.Sl_at, which represents a FER 
(FERONIA) protein kinase-like gene, and Cit.10594.1. 
Sl_at, which represents an EP3-like chitinase gene. 
Cit.ll529.1.Sl_s_at, Cit.11530.LSl, Cit.ll94.1.Sl_s_at, and 
Cit.3826.1.Sl_at were shown to be down-regulated by the 
Las infections in two reports [5,11] (Additional file 1). 
Another example is Cit.l923.1.Sl_s_at, which represents 
a protein degradation component similar to Arabidopsis 
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Figure 3 HLB disease response subnetwork. (A) The HLB defense subnetwork is displayed with three categories of major biological processes 
highlighted by different colors. (B) The close-up of the box in the subnetwork shown in Figure 3A. 
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CSN5A (COP9 Signalosome 5A) and is assigned the 
GO terms of both auxin response and defense response 
[20,21]. This hub is interconnected to at least two large 
defense hubs, Cit.4216.1.Sl_s_at (a putative monooxy- 
genase similar to Arabidopsis M02) and Cit.2848.1.Sl_at 
(a leucine-rich repeat family protein). However, there is 
one hormone (ABA) response hub, Cit.4553.1.Sl_s_at, 
which represents nicotinamidase similar to Arabidopsis 
NIC1 (Nicotinamidase 1) and is linked to at least three 
defense hubs. One of the connecting hub genes, 
Cit.23352.1.Sl_at, is closest to Arabidopsis RLP33 
(Receptor-like protein kinase 33), and another two 
hubs, Cit.l0594.1.Sl_at and Cit.21654.1.Sl_s_at, represent 
EP3-like chitinase genes. Although Cit.4553.1.Sl_s_at itself 
is not HLB responsive, the above three defense hubs to 
which it connects were reported to be up-regulated in 
some of the transcriptomic studies (Additional file 1). 

The finding that the defense and hormone hubs are 
intertwined or overlapped (Figure 3) indicates a potentially 
important role for hormones in the HLB response in citrus. 
Also given by the increasingly clear roles for some hor- 
mones such as ethylene, ABA, JA and SA in plant defense 
response, we decided to analyze in more detail the hor- 
mone response subnetwork. In the HLB response network, 
GO terms for the response to auxin, GA, ABA, ethylene, 
JA and SA are overrepresented (Additional file 5) based on 
the hypergeometric method provided in the agriGO web 
tool [22] and thus the nodes for these GO terms are color- 
coded in the hormone response subnetwork (Figure 4) and 
listed in Additional file 9. It should be noted that four 
(auxin, GA, JA and SA) of these six overrepresented hor- 
mone GO terms are also determined to be overrepresented 
by using several algorithms implemented in the R package 
topGO which are proposed to eliminate local dependencies 
between GO terms [23,24]. It has been demonstrated that 
SA signaling is important for both local disease resistance 
and systemic acquired resistance [25] and a recent report 
showed the success in engineering the NPRl-mediated SA 
signaling pathway to improve citrus resistance to another 
destructive disease canker [26]. Therefore, we used the SA 
response subnetwork as an example of performing the spe- 
cific hormone response network analysis. Using 49 SA re- 
sponse Probesets as the seed nodes, we constructed the SA 
response subnetwork consisting of 476 Probesets and 631 
interactions (Figure 5, Additional file 9). In the SA re- 
sponse subnetwork, there are two major subsets, each with 
several large hubs. The first major subset contains tran- 
scription factors similar to Arabidopsis AS1 (Cit.11529.1. 
Sl_s_at and Cit.ll530.1.Sl_at) and WRKY40 (Cit.10816.1. 
Sl_at), protein degradation component UBQ10 (Cit.107.1. 
Sl_x_at) and carbohydrate metabolic enzyme GSTU7 
(Cit.6954.1.Sl_at and Cit.ll419.1.Sl_at). The second major 
subset of the SA response subnetwork has two large hubs, 
both of which represent the UBQIO-like protein 



degradation component. A further analysis on this subset 
revealed that besides the two UBQ10 hubs, two other tran- 
scription factors closest to AS1 (Cit.42831.Sl_at) and 
MYB16 (Cit.l9398.1.Sl_at) serve as smaller hubs linking 
the larger UBQ10 hubs (Figure 6). WRKY, MYB and AS1- 
like transcription factors have been reported to play im- 
portant roles in Arabidopsis defense responses 
[18,19,27,28]. Ubiqutin-mediated proteasome has also been 
shown to be critical for plant disease resistance [29,30]. Ac- 
cumulating evidence suggest that WRKY, MYB and AS1- 
controlled transcriptional events and ubiqutin-mediated 
proteasomal degradation are critical for SA signaling 
[18,25,31,32]. Therefore, these results strongly indicate that 
protein degradation and transcriptional regulation are 
likely critical regulatory processes involved in SA response 
in citrus when challenged with the Las infection. 

Analysis of the early stage HLB response subnetwork 

At early stage, the HLB bacterium could rarely be 
detected, nor any HLB symptom observed, but the re- 
sponse to HLB in citrus could occur early at least at 
the transcriptional level. Therefore, we decided to 
analyze the subnetwork for the early-stage HLB re- 
sponsive genes. A total of 222 Probesets, including 158 
up-regulated and 62 down-regulated Probesets (Table 1, 
Additional file 10), were used as the seed nodes to 
map the HLB response network, resulting in the HLB 
early response subnetwork (Figure 7). This subnetwork 
based on the first degree neighbors of these seed 
nodes contains 461 Probesets and 683 interactions 
(Additional file 10). Among those Probesets, 29 are 
involved in carbohydrate metabolic process, 23 in nitrogen 
and amino acid metabolic process, 67 in transport, 27 in 
defense response, 24 in signaling and 24 in hormone re- 
sponse. GO enrichment analysis shows that carbohydrate 
metabolic process (p-value = 3.6E-4), transport (p-value = 
4.1E-10) and defense (p-value = 6.3E-3) are overrepre- 
sented. Although the hormone response category is not 
overrepresented, JA response consisting 10 Probesets is 
overrepresented with a p-value of 0.01). Therefore, our 
analysis of the early stage subnetwork indicates that even 
at this stage (5-9 weeks after inoculation), several im- 
portant biological processes have been activated or inacti- 
vated. In the HLB early response subnetwork, there is 
only one subset that has several large hubs, while all 
other small subsets have interactions that are not con- 
nected further (Figure 7). To provide further detail of the 
early stage response in citrus, we analyzed the two nodes 
in the large subset of this subnetwork, Cit.29252.1.Sl_s_at, 
and Cit.l2214.1.Sl_s_at. 

Cit.29252.1.Sl_s_at represents a triacylglycerol lipase 
gene most closely related to Arabidopsis EDS1 (Enhanced 
disease susceptibility 1) [33]. Extracting this EDSl-\ike gene 
from the HLB early response subnetwork (Figure 7) shows 
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Figure 4 HLB hormore subnetwork. The HLB hormone response subnetwork is displayed, with six hormone response groups coded with 
distinct colors. 

^ ) 



that EDS1 interacts with 15 Probesets (Figure 8A). Among 
these Probesets, one Probeset (Cit.38351.1.Sl_at represent- 
ing a MATE efflux family protein) has interactions with 
only two other Probesets, five Probesets (Cit.6535.1.Sl_at, 
Cit.l5432.1.Sl_at, Cit.l0234.1.Sl_s_at, Cit.4135.1.Sl_s_at, 
and Cit.2933.1.Sl_s_at) form the large hubs each with 50- 
113 interactions, and nine other Probesets form the 
medium-size hubs with 11-44 interactions. The fact that 
Cit.29252.1.Sl_at connects with the five large hubs indi- 
cates a potentially critical role in citrus response to the 



HLB bacterial infection. Cit.6535.1.Sl_at represents a 
carbohydrate transmembrane transporter or phosphate 
transmembrane transporter, Cit.l0234.1.Sl_s_at is closely 
related to CB5-E (Cytochrome B5 isoform E) involved in 
heme binding, Cit.4135.1.Sl_s_at represents a putative 
CC-NBS-LRR class disease resistance protein, and 
Cit.2933.1.Sl_s_at is very similar to Arabidopsis HMGB1 
(High mobility group 1) involved in transcriptional control 
through chromatin remodeling. In addition, some of the 
medium-size hubs that interact with the EDS1 -like gene 
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Figure 5 SA response subnetwor. The SA hormone response subnetwork is displayed. Nine large hubs are labeled with Probeset ID. 



play important roles in protein modifications or lipid me- 
tabolism. For example, Cit.39054.1.Sl_s_at is closely 
related to Arabidopsis SAG101 (Senescence-associated gene 
101) which encodes a carboxylesterase or triacylglycerol 
lipase and has been shown to physically interact with EDS1 
[34], and Cit.373.1.Sl_x_at is similar to UBQ10 (POLYU 
BIQUITIN 10). The EDS1 -like citrus gene was up- 
regulated at the early stage (2.0-fold difference) and at the 
very late stage (6.5-fold difference; Additional file 1) in only 
one of the four studies [11], and most of 15 hub genes that 
interact with the citrus EDS1 -like gene were also up- 
regulated by the Las infection in some of the studies with 
the exception of Cit.373.1.Sl_x_at (UBQ10), Cit.39054.1. 



Sl_s_at (SAG101) and Cit.l0182.1.Sl_s_at (an aldo/keto 
reductase family protein) (see Additional file 1). Therefore, 
the finding that so many HLB responsive hub genes in cit- 
rus connect to EDS1, which is critical for disease resistance 
in Arabidopsis and other plants [33,35,36], indicates that 
EDSl-mediated defense response mechanism might be im- 
portant in citrus response to the HLB bacterial infection at 
least at early stage. 

Cit.l2214.1.Sl_s_at represents a transcription factor most 
closely related to Arabidopsis NAC096. Mapping this Probe- 
set as the seed node to the HLB response network with the 
second degree neighbors resulted in an NAC096 subnet- 
work (Figure 8B). Two medium-size hubs were identified in 
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Figure 6 A subset of SA response subnetwork. A subset of the SA hormone response subnetwork with four large hubs is displayed. 



this subnetwork: Cit.l0032.1.Sl_x_at and Cit.l5606.1.Sl_at, 
both of which were up-regulated transcriptionally by the Las 
infection (Additional file 1). Cit.l0032.1.Sl_x_at represents a 
GA-responsive GAST1 homolog and has been reported to 
be responsive to other hormones such as BR and ABA 
[37,38]. Cit.l5606.1.Sl_at has interactions with 15 Probe- 
sets and is closely related to Atlg80130 which encodes 
Arabidopsis tetratricopeptide repeat (TPR)-like superfamily 
protein and is responsive to oxidative stress [39]. Given 
that both GA response and oxidative stress response have 
been implicated an important role in a relatively resistant 
variety US-897 in response to the Las infection at the 
very late stage [12], our preliminary analysis of the 
NAC096 subnetwork supports that transcriptional control 
involving hormone response and oxidative stress response 
might also be important even at the early stage of the 
HLB bacterial attack. 



Subnetwork analysis reveals transport process as a key 
component in the HLB response core subnetwork 

It is likely that the commonly up-regulated genes can de- 
fine a default or core response pathway for citrus plants 
to resist the attack by the HLB bacteria; we therefore 
attempted to address whether there is a common subnet- 
work that could be affected by HLB. We mapped 21 
commonly up-regulated Probesets into the HLB response 
network, resulting in the formation of the HLB core sub- 
network (Figure 9). This subnetwork based on the first 
degree neighbors contains 123 Probesets and 181 inter- 
actions (Additional file 11). The hub gene analysis shows 
that the subnetwork has eight large hubs, all of which 
were up-regulated, and four small hubs (Figure 9). 
Among six categories of biological processes analyzed in 
the HLB response network, transport and carbohydrate 
metabolic process were overpresented in this core 
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subnetwork (p-value = 1.0E-05 and 1.3E-4, respectively). 
The Probesets belonging to the three categories, carbo- 
hydrate metabolic process, transport and hormone re- 
sponse, were then color-coded in this subnetwork 
(Figure 9). The largest category of the genes in the sub- 
network belongs to transport process, with a total of 22 
Probesets (representing 18% of the nodes in the core sub- 
network versus 6% in the HLB response network). 
Among these Probesets, four form the large hubs: 
Cit.ll459.1.Sl_s_at, Cit.ll460.1.Sl_at, Cit.3171.1. 
Sl_x_at, and Cit.l7561.1.Sl_s_at. Given the importance 
of hub genes in the biological networks and overrepre- 
sentation of transport in the subnetwork, we propose 
that transport process is a key component in the HLB re- 
sponse core subnetwork. 

There are 13 Probesets grouped into the category of 
carbohydrate metabolic process and 11 Probesets that be- 
long to the hormone response category. For the category 



of carbohydrate metabolic process, Cit.l3437.1.Sl_s_at 
forms a larger hub with 11 interactions, and Cit.17155.1. 
Sl_at forms a smaller hub with seven interactions. 
Cit.l3437.1.Sl_s_at represents a citrus gene similar to 
Arabidopsis APL3 encoding a glucose-l-phosphate adeny- 
lyltransferase. Cit.l7155.1.Sl_at represents a gene closely 
related to BGLU11 (Beta-glucosidase 11) hydrolysis of O- 
glycosyl compounds. For the hormone response category, 
Cit.l9674.1.Sl_s_at forms a larger hub with 15 interac- 
tions, and Cit.l0032.1.Sl_x_at and Cit.25840.1.Sl_s_at 
form smaller hubs with seven and six interactions respect- 
ively. As described above, Cit.l9674.1.Sl_s_at represents a 
gene closely related to L0X2 encoding a lipoxygenase and 
exhibiting response to JA. In Arabidopsis, LOX2 has also 
been shown to be involved in JA biosynthesis in response 
to wounding [40,41] and recently in disease development 
[42]. As described previously, Cit.l0032.1.Sl_x_at repre- 
sents a GA-responsive GAST1 homolog and is connected 
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Figure 8 Early stage HLB response subnetworks based on two representative genes. Two gene-specific subnetworks are displayed with all Probesets 
labeled. (A) A subset of early stage HLB response subnetwork centered on the hub Cit.29252.1 .S1_s_at. Because many second degree interactions exist for 
most of the first degree nodes connecting to Cit.29252.1. S1_s_at, only the first degree interaction-based subnetwork is presented. (B) The NAC096 
transcription factor subnetwork is displayed. 
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Figure 9 HLB response core subnetwork. The HLB response core subnetwork is constructed using the 21 commonly up-regulated Probesets 
(see Table 2) as the seed nodes (highlighted in red). First degree interactions are shown. Green nodes indicate "transport", yellow nodes represent 
"hormone response", and blue nodes indicate "carbohydrate metabolic process". Nodes overlapped in the categories of "transport" and "hormone 
response" exhibit a mixed color of green and yellow. 



to the NAC096 transcription factor subnetwork in the 
HLB early response subnetwork (Figures 7 and 8B). Inter- 
estingly, Cit.25840.1.Sl_s_at represents a gene very similar 
to Arabidopsis WBC11 (White-brown complex homolo- 
gous protein 11) which encodes an ATPase coupled to 



transmembrane movement of substances or fatty acid 
transporter [43]. This small hub is responsive to ABA and 
salt stress but is also involved in fatty acid transport, imply- 
ing a potential role for hormone signaling in the control of 
transport process. 
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The remaining two large hubs in the HLB response 
core subnetwork are formed by Cit.l2172.1.Sl_s_at and 
Cit.l5630.1.Sl_at. Cit.l2172.1.Sl_s_at represents a puta- 
tive O-methyltransferase family 2 protein most closely 
related to the protein encoded by At4g35160. At4g35160 
is only annotated as a general GO term "methylation", 
and predicted to contain a winged helix-turn-helix tran- 
scription repressor DNA-binding domain without any 
functional implication. This hub includes 31 interac- 
tions, and most of the interactions are with the Probe- 
sets related to transport process. Cit.l5630.1.Sl_at 
represents a gene closest to At4g33040 which encodes a 
glutaredoxin family protein. It connects to a transportor 
hub through Cit. 17265. l.Sl_at (a glycine-rich protein 
without any specific functional annotation) and the two 
hormone response hubs through Cit.l7398.1.Sl_at 
(which represents e gene encoding an unknown protein). 
In Arabidopsis, At4g33040-encoded glutaredoxin family 
protein is annotated as a protein involved in cell redox 
homeostasis, which could have a potential role in re- 
sponse to oxidation stress. 

Analysis of a phloem protein subnetwork implicates a 
potential role for zinc transport in the citrus HLB defense 
response 

Given the potential importance of phloem protein 2 
(PP2)-type lectin in phloem morphogenesis in particu- 
lar the formation of sieve plug [5,6,44], PP2-like pro- 
tein genes in citrus were used as an example to 
further illustrate the application of the HLB response 
network. A survey of ten PP2-like genes present in the 
citrus GeneChip showed that four of the PP2-like 
genes were up-regulated and one down-regulated 
(Additional file 12). Although their expression pattern 
was quite different, one gene represented by the Pro- 
beset Cit.35955.1.Sl_at was dramatically up-regulated 
at late stage (44.6-fold difference) and very late stage 
(42-506 fold differences) in all of the four reports 
[5,10-12] except for the relatively resistant variety US- 
897 which did not exhibit any activation at very late 
stage (Additional file 1). This gene is closely related to 
Arabidopsis PP2-B8. This Probeset and the other, 
Cit.3272.1.Sl_s_a, are present in the HLB response net- 
work. The latter one represents a PP2-A15-like gene but 
expression of this gene was not affected by HLB in any of 
the four reports and it only connects with three genes in 
the HLB response network (Additional file 12). The lack 
of activation of Cit.35955.1.Sl_at by the Las infection at 
the early stage might be due to that the HLB symp- 
tom has not been fully developed yet. When the PP2- 
B8 subnetwork was constructed, we found that this 
gene connects with 20 Probesets (Figure 10A) which 
are interconnected frequently between each other 
(Figure 10B). Furthermore, seven of the 20 first degree 



interacting Probesets represent the genes involved in trans- 
port, and three of these genes are predicted to encode zinc 
transporters (Additional file 12). In addition, four Probesets 
represent genes encoding zinc-binding proteins. Given that 
HLB disease symptom was initially thought to be related to 
zinc deficiency [1,2], our network analysis approach pro- 
vides an intriguing possibility for zinc transporters or zinc- 
binding proteins to function in citrus response to the HLB 
bacterial infection. 

Discussion 

The transcriptomes in citrus in response to the HLB 
bacterial infection have been well documented in four 
previous reports [5,10-12], but the information regarding 
the interactions between the differentially expressed 
genes is lacking. Through the combination of transcrip- 
tome comparative study and gene coexpression network 
analysis, we have provided for the first time a systems 
view of how the citrus host plant exerts a genome-wide 
response to the HLB bacterial infection. 

First, we have constructed an HLB response network 
involving 3,507 Probesets with 56,287 interactions. Using 
the transcriptome datasets and orthology-based or ex- 
perimentally verified protein-protein interaction data- 
sets, gene-gene interactions or interactomes have been 
constructed in the model plants including Arabidopsis 
and rice (for example, [16,45-52]) and occasionally in 
non-model plants such as soybean [53] and barley [54]. 
However, there was no report on gene-gene interaction 
networks in citrus prior to our work. We used the Pec 
method to construct a gene coexpression network in cit- 
rus, with a focus on the HLB response mechanism. The 
citrus gene coexpression network will be very useful for 
the citrus researchers to visualize the subnetworks spe- 
cific for certain biological processes (such as carbohy- 
drate metabolic process or signaling not analyzed in 
detail in this report), or to search some potential gene- 
gene interactions for certain genes or a group of genes 
in the future. The Citrus Gene Interaction Networks 
(CitGIN) database has been constructed and made avail- 
able to the research community to query through the 
Internet (http://xt.cric.cn/ cgr/ CitGIN.php) . 

Second, our analysis of the defense subnetwork has 
shown that many defense hubs and hormone hubs are 
intertwined or overlapped. Although the roles of hormone 
and defense response genes have been discussed in the 
four previous reports [5,10-12], our network analysis fur- 
ther indicates that hormone response is interconnected to 
defense response in citrus when challenged by the HLB 
bacteria. This may lead to the development of integrating 
hormone and disease response pathways as a potentially 
more effective genetic means to improve the citrus resist- 
ance to HLB. 
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Cit.34982.1.S1_at Cit.3611.1.S1_s_at 




(PP2-B8) 

Figure 10 Phloem protein (PP2) subnetwork. The subnetworks for the PP2, based on first (A) and second (B) degree neighbors of Cit.35955.1. 
S1_at. Six of the Probesets belonging to the second degree interactions (B) are labeled with yellow. 
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Third, our comparative studies of transcriptomes have 
led to the identification of subsets of commonly up- 
regulated and stage-specific HLB responsive genes. In 
contrast to those four GeneChip reports where various 
statistical methods and fold-change cutoffs were used, 
we used the same procedure for the analysis of all of the 
transcriptome datasets. Furthermore, by mapping the 
subset of commonly up-regulated genes into the HLB 
response network, we have found that the genes belong- 
ing to the categories of carbohydrate metabolic process, 
transport and hormone response are positioned as the 
large hubs in the HLB response core subnetwork. This 
indicates that these three processes constitute a core 
subnetwork for the citrus host response to the HLB bac- 
terial infection. In addition, we propose that transport is 
a key component in this HLB response core subnetwork. 

Fourth, using PP2 gene as an example of applying the 
HLB response network, our subnetwork analysis pro- 
vides an intriguing possibility for the zinc transporter or 
zinc binding proteins to act with PP2 protein in re- 
sponse to the HLB bacterial infection. PP2 proteins be- 
long to a large gene family in higher plants. However, 
they have not been assigned a specific biological process, 
and thus their biological function remains unknown. 
They are predicted to bind carbohydrates and have been 
implicated a role in the formation of sieve plug or re- 
placement phloem [5,10,12,44]. Some of the PP2 genes 
from other organisms such as melon, cucumber and 
Arabidopsis are specifically or preferentially expressed in 
companion cells but their protein products are translo- 
cated in sieve elements [44]. This indicates a role for 
PP2 proteins not only in intracellular signaling but also 
in long distance intercellular communication [44]. Re- 
cent evidence show that PP2-type proteins might be 
involved in aphid-mediated virus transmission in cu- 
cumber [55] and overexpression of PP2-A1 in Arabidop- 
sis increased resistance to a phloem-feeding insect [56], 
indicating a possibly more active role for PP2 family pro- 
teins in defense response in plants. In citrus, the Probe- 
set Cit.35955.1.Sl_at, which is closely related to 
Arabidopsis PP2-B8, was dramatically up-regulated at 
late stage and very late stages. The most surprising fea- 
ture of the PP2-B8 subnetwork is that the 20 Probesets, 
which are the first degree neighbors of Cit.35955.1. 
Sl_at, are interconnected frequently between each other. 
This indicates that these genes might be regulated by the 
precise coordination of various signaling pathways 
through transcription factors, chromatin modification or 
remodeling proteins or other factors. Furthermore, seven 
of the 20 interacting Probesets encode proteins involved 
in transport, consistent with our proposal that transport 
is a key component in the HLB response core subnet- 
work. In addition, three of the seven transporters are 
predicted to transport zinc, and the PP2 subnetwork also 



contains four Probesets which represent the genes en- 
coding zinc-binding proteins. Intriguingly, HLB disease 
symptom was initially thought to be related to zinc de- 
ficiency and the zinc transport system is required for 
virulence in other organisms [1,2], and therefore the 
PP2 subnetwork analysis indicates that zinc transpor- 
ters or zinc binding proteins may have a potentially 
important role for citrus to respond to the HLB bac- 
terial infection. Taken together, our analysis using the 
HLB response network can lead to an intriguing but 
testable hypothesis (which cannot be done without the 
network analysis) regarding the role of PP2 proteins 
and zinc transport system or zinc binding proteins in 
citrus HLB defense response. 

It should be noted that there are some potential limita- 
tions in our network study. The first one is GO enrich- 
ment analysis. The agriGO web tool, which is based on 
the hypergeometric method and used in this work, does 
not take into account the local dependency of GO terms. 
Using the four algorithms provided in the topGO R pack- 
age which are proposed to eliminate local dependencies, 
we have found that four of the six hormone GO terms 
determined to be overrepresented by agriGO are also 
overrepresented, while the two other hormones (ABA 
and ethylene) have their child GO terms being truly over- 
represented. Therefore, different algorithms or statistical 
methods in GO enrichment analysis will probably lead to 
some differences in terms of the overrepresented GO 
terms for the nodes in the HLB response network. 

The second limitation is due to the small sample size. 
Computational prediction of gene-gene interactions usu- 
ally requires large sample size; however relatively small 
number of samples were recently used to construct gene 
coexpression networks specific to certain aspects of biol- 
ogy (such as [50,52,57]). In our analysis, we used the 
transcriptome datasets described in four previous 
reports [5,10-12]. Among these, only one study is avail- 
able for early stage and late stage respectively, while 
there are four studies for very late stage (Table 1). We 
have found that only a small number (222) of signifi- 
cantly regulated Probesets can be identified for early 
stage, while almost 600 and 2,000-4,000 differentially 
expressed Probesets can be found for late and very late 
stages respectively. The variation in the number of dif- 
ferentially expressed genes at different stages could be 
caused by the difference in experimental conditions 
given that different ages and varieties of trees and differ- 
ent sources of inoculants were used in different years in 
those four reports. However, this variation might lead to 
some sort of bias towards the very late stage genes. To 
minimize the possibility that the interactions we have 
detected were the result of random events due to the 
small sample size, we have selected a high Pec cutoff 
value which has led us to believe that the interactions 
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are more likely statistically significant rather than by 
random and that the topology of the HLB response net- 
work is quite similar to most biological networks. Fur- 
thermore, the cross-validation result shows a high 
degree of preservation of gene coexpression patterns, 
suggesting that the HLB response network is at least 
moderately robust and biologically relevant. Therefore, 
despite some limitations due to the small sample size 
and the experimental variations, the network reported 
here should be quite useful for the citrus research com- 
munity and have provided some novel insights into the 
citrus HLB defense mechanisms. When larger scale tran- 
scriptome datasets become available in the future, simi- 
lar network analysis will provide a comprehensive 
picture of the gene networks in citrus. 

The most daunting challenge in the citrus post- 
genomic era remains how to identify the best candidate 
genes for functional dissection of the HLB response 
mechanism and for genetic modification with an ul- 
timate goal of improving the HLB resistance in 
citrus. Genetic variations of HLB susceptibility [9] 
clearly shows the potential towards dissection of gen- 
etic mechanisms of HLB resistance, but understand- 
ing the inheritance patterns and subsequently cloning 
the disease genes requires a long term effort because 
of long juvenile phase and complex reproductive 
biology for citrus. Recent developments have shed 
some lights into the identification of key hub genes 
as candidate regulatory genes. For example, a seed 
germination study found that 22-50% of the Arabi- 
dopsis hub genes identified from the seed germin- 
ation network actually have physiological functions in 
the control of seed germination [49]. Therefore, the 
hub genes identified in this report may potentially be 
the first batch of candidates for the functional test in 
HLB resistance in citrus. 

Conclusions 

Through integration of transcriptome comparison and 
gene coexpression network analysis, we have provided 
novel insights into the mechanism by which citrus host 
plants respond to the HLB bacterial infection. Specifically, 
several biological processes are important in the citrus 
HLB response network, including carbohydrate metabolic 
process, nitrogen and amino acid metabolic process, 
transport, defense response, signaling and hormone re- 
sponse. Furthermore, our results have led us to propose 
that transport is a key component in the HLB response 
core subnetwork. This systems view of citrus response to 
the Ca. Liberibacter spp. infection will be a critical first 
step towards dissecting the genetic mechanisms of HLB 
response and ultimately improving HLB resistance in 
citrus. 



Methods 

Data collection and preprocessing 

Raw data for citrus Affymetrix GeneChip analysis pub- 
lished by Fan et al. [11] (GSE29633) and Albrecht and 
Bowman [12] (GSE30502) were downloaded from NCBI. 
Raw data published in [10] (GSE33459) and [5] were 
kindly provided by Drs. Bowman and Wang, respectively. 
These .eel files were read into R and preprocessed using 
"rma" function [58] and normalized using the "normalize, 
quantiles.robust" function. After quantile normalization, 
Probesets with an absent (A) call were removed using the 
"pma" function. Probesets with the calls of present (P) or 
marginal (M) in at least two samples in each of the four 
reports above were included in the analysis. All of the stat- 
istical analysis and gene expression network construction 
were performed in the R environment. 

Analysis of significantly regulated genes 

The adjusted local pooled error (LPE) method [13] was 
used to identify differentially expressed transcripts, as this 
method has been shown to provide high power in analyz- 
ing microarray data with small sample size. A gene was 
called statistically significant if its permutation-based false 
discovery rate (FDR) p-value was smaller than 0.05 and at 
least a two-fold change was observed. 

Network construction and visualization 

For computational reasons, up to 10,000 of the Pro- 
besets with highest expression levels were selected 
from each of the datasets described in the four 
reports [5,10-12]. The HLB responsive genes identi- 
fied in this study (described above) were then added 
to this list and duplicated ones were removed, result- 
ing in a total of 10,668 common Probesets for each 
of the four datasets. Gene coexpression network was 
constructed from the preprocessed files using R 
package "weighted correlation network analysis" [59]. 
Following the protocol for constructing gene co- 
expression network using multiple datasets [15], we 
first calculated Pearson correlation matrix for each 
dataset. We then obtained an overall weighted correl- 
ation matrix based on the number of samples used in 
that dataset. The weight for each correlation matrix 
was defined as Wi = . — \^HVw where «,• was the 

V" Hi + 5 V" max 

i 

number of samples for ith dataset, n max was the maximum 
number of samples in all datasets, and s was the number 
of datasets used. Two nodes were determined to be con- 
nected if the absolute value of the Pearson correlation 
coefficient (Pec) exceeded 0.93. The threshold of 0.93 
was selected such that it gave the best overall fit to each 
dataset based on the criteria such as the scale-free top- 
ology model fitting index, mean network connectivity, and 
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network density. Cytoscape v2.8.2 [60] was used to 
visualize the networks and Photoshop was used to edit the 
images. 

GO analysis and Arabidopsis orthology prediction 

Because of the lack of citrus genome annotation for 
the Probesets in the Affymetrix chip, the Probesets 
were used for all analysis. They were annotated using 
Arabidopsis orthologs or homologs. The Probesets 
were annotated by searching against the Arabidopsis 
genome using the tool provided in HarvEST database 
(http://harvest-web.org). GO terms were assigned to 
the citrus Probesets based on their corresponding Ara- 
bidopsis gene ID (AtGID). For those without AtGID, 
general GO terms were assigned: "biological process" 
(GO:0008150), "molecular function" (GO:0003674), 
and "cellular component" (GO:0005575). GO enrich- 
ment analysis was performed using the hypergeometric 
statistical method with Hochberg FDR adjustment in 
the AgriCO website (http://bioinfo.cau.edu.cn/agriGO/ 
analysis.php) as described elsewhere [22]. 
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